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We present a theory for the hnear dynamics of a weakly interacting Bose gas confined inside a 
harmonic trap at finite temperature. The theory treats the motions of the condensate and of the 
non-condensate on an equal footing within a generalized random-phase approximation, which (?) 
extends the second-order Beliaev-Popov approach by allowing for the dynamical coupling between 
fluctuations in the thermal cloud, and (ii) reduces to an earlier random-phase scheme when the 
anomalous density fluctuations are omitted. Numerical calculations of the low-lying spectra in the 
case of isotropic confinement show that the present theory obeys with high accuracy the generalized 
Kohn theorem for the dipolar excitations and demonstrate that combined normal and anomalous 
density fluctuations play an important role in the monopolar excitations of the condensate. Mean- 
field theory is instead found to yield accurate results for the quadrupolar modes of the condensate. 
Although the restriction to spherical confinement prevents quantitative comparisons with measured 
spectra, it appears that the non-mean field effects that we examine may be relevant to explain the 
features exhibited by the breathing mode as a function of temperature in the experiments carried 
out at JILA on a gas of *^Rb atoms. 

PACS numbers: 03.75.Kk, 05.30.Jp, 67.40.Db 

I. INTRODUCTION 

Soon after the realization of Bose-Einstein condensation in trapped atomic gases, an important development in 
this field has been the measurement of the frequencies and damping rates of collective excitations Q, l^i Hi ^ IS- 
These measurements are very accurate and provide a unique opportunity for quantitative tests of quantum theories 
of the dynamics of many-body systems. In particular, the measurements of the lowest-energy excitations made at 
JILA 0] on ^''Rb gases at various temperatures have proved hard to understand at simple mean-field level 0, Q 
and have therefore stimulated a number of theoretical studies to address effects beyond the mean-field approximation 
MIIiLl, 12, 13, ,14, 15, 16]. 

The key issue in investigations transcending the mean-field level is the full dynamic description of both condensed 
and non-condensed atoms and their mutual interactions 0. While the condensate dynamics is well described by 
a single nonlinear Gross-Pitaevskii equation (GPE), how to monitor the evolution of the non-condensate is a much 
more delicate problem. The best candidate theory that takes into account the coupled dynamics of condensate and 
non-condensate for a homogeneous weakly-interacting Bose gas in the collisionless limit is the second-order Beliaev- 
Popov (SOBP) theory jT^, which has been reexamined recently by Shi and Griffi n 1 181 and extended to trapped 
gases by Fedichev and Shlyapnikov and by Giorgini T?] (see also Rusch et al. jl5|). However, for the trapped 
gas the Thomas- Fermi approximation on the SOBP theory fails to account for the JILA observations [loL H^. One 
possible reason is that the dynamics of the condensate and non-condensate are not treated on an equal footing in 
the theory, i.e. the dynamical coupling between fluctuations in the thermal cloud is not included. This coupling 
should be important when the thermal fraction is significantly populated and, as will be discussed below, is in fact 
needed to satisfy the generalized Kohn theorem for the dipole modes. One way to include these processes is to use 
the linear response theory in the random-phase approximation (RPA) as developed by two of us ■ Such a treatment 
chooses the Hartree-Fock gas as the reference system for the thermal atoms, thus neglecting the anomalous density 
fiuctuations that may play a role at intermediate temperatures. 

In the present paper we improve on the Hartree-Fock RPA (HF-RPA) by including the anomalous density fluctua- 
tions. The resulting theory can be referred to as the HFB-RPA since our choice of the reference system is provided 
by the first-order Hartree-Fock-Bogoliubov theory. We explicitly show that the HFB-RPA theory formally reduces 
to the SOBP theory given by Fedichev and Shlyapnikov ^3 ^-nd by Giorgini if (i) one excludes the process of 
driving the non-condensate by its self-generated dynamical potential, and (ii) one keeps only terms of second order in 
the coupling constant. It is interesting to note that the HF-RPA similarly reduces to the dielectric formalism given 
by Reidl et al. jl^. 
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We then numerically investigate the low-lying excitations of a fluid representing a Bose-condensed gas of 2000 
^^Rb atoms in a spherically symmetric harmonic trap at finite temperature by using the HFB-RPA as well as the 
SOBP theory and the HF-RPA. All three theories give qualitatively the same results for the quadrupolar mode of 
the condensate. However they predict different trends for the monopolar mode, due to the strong coupling between 
the oscillations of the condensate and those of the non-condensate. For the first time this observation highlights the 
crucial roles played already in the linear excitation spectra by the normal and anomalous density fluctuations of the 
non-condensate. 

The paper is organized as follows. In Sec. II we derive the generalized RPA equations within the framework of 
the Hartree-Fock-Bogoliubov approximation, and in Sec. Ill we briefly demonstrate how to deduce the SOBP theory 
from the HFB-RPA equations. In Sec. IV we describe our numerical procedure for calculating the spectral response 
functions and check their accuracy, and in Sees. V and VI we present our numerical results for the low-energy 
excitations. Finally, Sec. VII presents our main conclusions. 



II. THE HFB-RPA THEORY 



The essential idea of the RPA is that the gas responds as a reference gas to self-consistent dynamical potentials Q . 
In the HF-RPA treatment one chooses as dynamical variables the density fluctuations Sric of the condensate and 6n of 
the non-condensate The HF-RPA equations follow by imposing that the condensed and non-condensed particles 
experience dynamical Hartree-Fock potentials generated by both types of density fluctuations and respond to them 
as a Hartree-Fock gas. 

Our starting point for the derivation of the HFB-RPA is the definition of the appropriate single-particle reference 
system. The contribution of the anomalous density is included by choosing the Hartree-Fock-Bogoliubov gas at finite 
temperature as reference, which is defined in terms of the condensate wavefunction $o and of the single-particle 
amplitudes uj and Vj for the non-condensate The condensate is described by the generalized GPE 
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where we adopt the standard contact-pseudopotential model characterized by the coupling constant g — ATrh^a/m, 
with a being the s-wave scattering length. In Eq. Vext (r) = m{uj'^x'^ +uj'^y^ -|-ci;^z^)/2 is the external confinement 
2 ^o^.^_Y- , \„,.l^^\Af. , l,.^.m21 and to" (r) = V . [(1 + 2/,) (r) «* (r)] are 
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the condensate density and the normal and anomalous thermal densities, fj — 1/ (e'^^J — l) being the Bose-Einstein 
distribution with /3 = l/ksT and /i the chemical potential. The non-condensate amplitudes are obtained by the 
solution of the generalized Bogoliubov-deGennes equations 



C (r) u, (r) + g ($2 (r) + (r)) v, (r) = e,u, (r) 

C (r) Vj (r) + g (r) + rfi°* (r)) (r) ^ -e^v^ (r) . 



(2) 



Here C (r) = —Ti'\/^l1m -I- V^xt (r) + 2^ (ric (r) + (r)). The Popov approximation to the Hartree-Fock-Bogoliubov 
(HFB-Popov) theory is recovered by setting (r) = in Eqs. Q and ^ (l9l |. 

We would like to remark that from a dynamical point of view the amplitudes Uj and Vj can alternatively be viewed 
as excitations out of the condensate. The duality of such mean-field description follows from the assumption of Bose 
symmetry breaking (see e.g. 20] for a discussion). 

In deriving next the HFB-RPA equations we adopt five dynamic variables, which are the fiuctuations and 
(5$* of the condensate wavefunction and its complex conjugate, the normal density fluctuation Jn, and the anomalous 
density fluctuation 8m together with its complex conjugate brh* . 5$ and b^* are separately introduced because of their 
different coupling to bfn and (5m*, and are related to the density fluctuation of the condensate by bric = $q(5$-|-$o'5$*. 
The HFB-RPA then follows naturally by evaluating the self-consistent dynamical Hartree-Fock-Bogoliubov potential 
generated by the density fluctuations of the condensate (phonon quasiparticles) and of the non-condensate (thermal 
quasiparticles). This can be done by invoking the decomposition 



V'(r,t) = $0 (r)+V^(r,i) 
for the Bose field operator in the interaction Hamiltonian, 
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(3) 



H„t = I / dri^+(r,i)V'+(r,t)i^(r,i)^(r,i) 
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$o|^ + 2|$o|'$SV^ + 2|$o|'$oV3+ 
+ 4 |$o|' + $o$o^+^+ 



(4) 



Note that in the choice made in Eq. (jS}, which is different from those generally used in the literature, the non- 

equilibrium statistical average (^ijj{r,t)^ of the operator '0(r,t) is non-zero since we prefer to extract from rp(r,t) 

a time-independent condensate wavefunction. Rather, ^ (r, t) gives the field operator for the phonon quasiparticles 

and describes the condensate fluctuation, (^-tp (r, t)^ = {ip (r, t)) — $o (r) ~ $ (r, t) — $o (r) — <5$ (r, t) . Analogously 

,5$*(r,t) = (7A+(r,t))-$5(r). 

The sclf-consistcnt dynamical potentials are originated from the higher- order correlation terms beyond the mean- 
field description and are contained in the last line of Eq. |0J. We approximate these terms by using Wick's theorem 
in the following manner: 



4$;^ (ij+^pj ip + 2% l^^iP^ iP+ + A(^15^^P+^ + 2<^l5<^*^iP, 
4<I>o (v^+^a\ + 2$o ('^+^+\ + 4$o'5$*V'^V' + 2^o5^'^P+'^P+ 



and 



■ip+-ip+-ip-ip ~ 4 l^i^^'ip^ 'ip'^'ip + ijtp + (^-ipip^ 'ip^'ip^ . 



(5) 
(6) 

(7) 



Explicitly, the fluctuations of the non-condensate are defined by (5n (r, i) — {ip^ {r , t) ip {r , t)) — n"(r), 6rh{r,t) — 
{ip (r, t) ip (r, t)) — rhP (r), and 6m* (r, t) = {ijA (r, t) (r, t)) — mP* (r). We insert these definitions into Eqs. 
^ , remove the terms that are proportional to (r) , rhP (r) and mP* (r) as these are already accounted by the 
Hartree-Fock-Bogoliubov mean-field equations, and finally collect together the remaining terms. We then find that 
the self-consistent dynamical potential induced by fiuctuations is 

,sc 



JV""^ = g dr 



2%Sh-ip + ^QSrhip+ -\- 2$o^^V'^ + ^oSrh*ip 

+2%S^iP+iP + $o(5$*?/iV' + 2^oS^*iJ^iJ + ^aS^i'^i'^ 
-\-26n-tp+-tp -\- 6m-tp+-tp+ /2 + Sm*'ip'ip/2 . 



(8) 



Physically, the eight leading terms on the right-hand side of Eq. ijHJl are the self-consistent potentials generated by 
phonon quasiparticles on themselves and on thermal quasiparticles. These terms have been discussed by Giorgini |l4| 
and by Liu and Hu and, as we shall see explicitly below, lead to the SOBP theory in a perturbative treatment to 
second order in the coupling constant. On the other hand, the last three terms in Eq. © describe the self-potential of 
the thermal quasiparticles and are expected to excite zero-sound-like collective modes of the non-condensate. Although 
these terms are only of third order in the coupling constant and therefore are missing in the SOBP theory, they may 
have a significant role when the depletion of the condensate is large. They are also required for consistency with the 
generalized Kohn theorem. 

With the self-consistent Hartree-Fock-Bogoliubov potential in Eq. © and using the notation xf ^ 
/ dr'x (r, r', 1^) /(r'): can write the coupled HFB-RPA equations for (5$, (5$*, Sh, 5rh and 5m* in a compact 
matrix form. Within the linear response framework we have 



5^ 
5^* 
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Xc 



2%5h + $o<5to* 
2<^o5h -I- $gJm 



and 




X"n Xnfh Xnm+ 
Xfhh Xmm Xmrh+ 
Xm+n Xm + m Xrh'^rh^ 



2$o'5$* + 2(5n 
5m* /2 
^ofJ* + 5m/ 2 



(9) 



(10) 



n, m, or m"*") are the two-particle response functions of the 



In these equations XajS (a, /3 = c or c ) and Xab (a, b 
condensate and non-condensate components, respectively. They can easily be evaluated by using the quasiparticle 
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amplitudes obtained from the Hartree-Fock-Bogoliubov solutions in the standard finite-temperature Green's functions 
technique 22] . For the condensate we have 



and 

Xc-c- (r, r , a;) _ 2^ 1, 7,^.+ - ?ic.+ + 6, j ' ^^^^ 

where w"*" = cj + i?7 with 77 = 0+. The expressions for the two-particle response functions of the non-condensate are 
lengthier and we list them in Appendix A. 

The coupled HFB-RPA equations and (|10|l are the central result of this work. They reduce to the HF-RPA 
equations if one omits the anomalous density fluctuations of thermal quasiparticles. That is, the HF-RPA gives 



5nc{Y,Lu) = 2g J dr'xc (r, r'; w) Jn(r', w) (15) 

and 

Sn (r, Lu) = 2g dv'xnh (r, r'; uj) [5nc (r', w) + (5n(r', uj)] , (16) 



where Xc(r,r';u;) = %{v)xcc%{-r') -t- $0 (r) XEc$S(r') + *S (r) Xce*o(r') + $0 (r) XeE$o(i-')- One must accordingly 
take the Hartree-Fock reference system in the calculation of the response functions ^ . 



III. REDUCTION TO THE SECOND-ORDER BELIAEV-POPOV THEORY 

In this Section we show that the coupled HFB-RPA equations for the normal modes of the condensate simplify to 
those obtained in the SOBP theory if we neglect the self-coupling of density fluctuations in the non-condensate and 
keep only terms up to second order in the coupling constant g. This discussion also allows us to define a RPA form 
of the SOBP theory, that will later be used in our numerical calculations. 

If we neglect the terms in 5n, 6fh and 5rh* on the right-hand side of Eq. p(J|) and substitute this equation in Eq. 
0, we immediately obtain the self-consistent equations for the fluctuations of the condensate as 

^S<P \ _„2f Xcc Xc-c\^( 5'^ 

where the matrix T) is defined as 

Y Am+n Am+m Am+m+ J \ ^ y 

Equation (|17|l is already of second order in g and we shall regard it as providing a second-order Beliaev-Popov theory 
within a random-phase framework (SOBP-RPA). 

The SOBP-RPA differs only slightly from the SOBP theory presented in Ref. in the sense that it still keeps 
a class of terms beyond second order. In fact, to second order in the coupling constant we can describe the small 
oscillations of the condensate by a set (mqsc, Vqsc) of quasiparticle amplitudes with excitation energy Cosc- By setting 
((5$, (5$*) = {uosc^Vosc) in Eq. (|17|) and using Eqs. (|ll|) - (|14|l we find the eigenfrequency of the oscillations of the 
condensate as eosc + — 17, where 

5E-i^ = g" ( vl,,, ) V ( 
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g / dr$o [2 + Ksc) + Ksc^m + v*,Jm*] . (19) 
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In recent work two of us |2l| have explicitly shown that Eq. (|19|l agrees with the result for the eigenfrequency shift 
given by the SOBP theory of Giorgini ,14 |. 

IV. NUMERICAL PROCEDURE 

We turn to numerical illustrations of the excitation spectra with the main aim of comparatively examining the three 
theoretical approaches that we have introduced in Sees. II and III. We do this in the case of a spherically symmetric 
trap in view of the complexity of the calculations involved. 

We excite density fluctuations by applying a time-dependent perturbation of the form 
F{t) oc exp (iujt) J drVp (r) i/;'^ (r) ij; (r) . In the HFB-RPA this corresponds to adding the terms 

[XccVp^ + XccVp^o, XccVp^ + XccVp^o] and [xnnVp, Vp \ on the right-hand side of Eqs. 

© and (lion , respectively. The various density fluctuations are then calculated by the method of Capuzzi and 
Hernandez 23j . with a discretization of the dynamical equations on a spatial mesh of up to 256 points. The 
frequencies of the collective excitations of the system can be extracted from the resonances of the spectral function 
X (w), which is also the quantity of experimental interest. This is defined in the HFB-RPA as 

X {^) = Xc (^) + Xt i^) (20) 

where 

Xc H - -^Im J drVpir) ($*<5$ + $o'5$*) (21) 

and 

x'-^ (uj) = -i/m J drVp{r) {Sh + 6m + 6m*) . (22) 

Here the indices C and T refer to the contributions from the condensate and from the non-condensate. Other quantities 
of interest are the density fluctuations of the condensate and the non-condensate, which are readily extracted from 
the solution of Eqs. ® and ((Tn|l . 

The main technical difficulty in the numerical calculations is how to renormalize the ultraviolet divergence caused 
by the use of contact interactions p^ . The divergence appears in the equilibrium anomalous density to" (r) and 
in the response functions Xmm,+ ^-nd Xm+m- The simplest way to implement renormalization is by removing the 
zero-temperature component of the above quantities. This procedure is not fully correct as it neglects the quantum 
contributions but these are extremely small at temperatures where the thermal corrections become important. 
Alternatively one can apply renormalization by regularizing to" (r), Xrhm+ Xm+m in real space |25|. We have 
checked that these two procedures give almost the same mode frequencies in calculations based on the standard SOBP 
theory. 

In brief, the numerical method that we have used consists of three steps. Firstly, we solve the HFB Eqs. ^ and 
(121) (or, in case of HF-RPA, the corresponding HE equations) to determine the equilibrium densities and quasiparticle 
amplitudes. We then construct the bare two-particle response fimctions and compute the dynamic fluctuations from 
Eqs. and H10|l . We finally calculate the imaginary part of the response functions according to Eq. (|20|l . In the 
present case of an isotropic trap, the calculations can be greatly simplified by projecting the RPA equations and 
the response functions onto the various multipole modes |23|. We shall be interested in the monopolar, dipolar, and 
quadrupolar excitations, which require setting Vp (r) oc r^, Vp (r) oc rcosO and Vp (r) oc r^F20 (^7 '■p)- 

In the following we evaluate a gas of iV = 2000 ^^Rb atoms in a spherical trap with trap frequency ujq = 2n x 182.5 
Hz, this value being the geometric average of the axial and radial frequencies in the JILA experiments 0. The 
temperature is taken in units of the critical temperature for an ideal gas with the same value of N and wq, which is 
Tc = 0.94?ia;oA^"'^/^. In most calculations we use a basis of n < rimax = 24 and I < Imax = 32 for the quasiparticle 
wavefunctions, where the indices n and / label the number of radial nodes and the orbital angular momentum of the 
wavefunction. 



A. Tests of numerical accuracy 

In this subsection we report some tests of the accuracy of our numerical calculations. First of all, we must replace 
the positive infinitesimal quantity 77 in the reference response functions by a finite value. In Fig. 1 we show the spectral 
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FIG. 1: Spectral responses (in arbitrary units) for the monopolar excitation as functions of frequency u (in units of uoo) as 
calculated from the HF-RPA at T /Tc = 0.5, plotted for two values of r) (in units of luo) as indicated in the panels. The three 
panels display the total spectral response (a) and the contributions of the condensate (b) and of the non-condensate (c). 




FIG. 2: Spectral response (in arbitrary units) for the monopolar excitation as a function of frequency u (in units of tJo) as 
calculated from the HF-RPA at T/Tc — 0.5 and rj = 0.005ljo with two kinds of basis set. 



functions for the monopolar excitation in the HF-RPA for two values of 7y at a reduced temperature T/Tc — 0.5. For 
small value of rj many spikes appear in the spectrum, due to the discrete basis set that was chosen for the dynamical 
description. With increasing rj these spikes are rounded off into broad resonances, which are insensitive to the precise 
value of rj. In the following we preferentially take 77 = O.Olwo in calculating the spectral functions, this choice being 
consistent with a typical experimental energy resolution 0. 

The other aspect of the calculations that needs examining is the role of the basis set. In Fig. 2 we show the HF- 
RPA monopole spectrum at T/Tc — 0.6 and rj — O-OOSljo, as calculated from two choices of basis set. These are the 
standard set as described above (solid line) and a set in which the number of basis function has been doubled (dashed 
line). No quantitative changes are found for the condensate response around oj — 2.2ujo, while for the response of the 
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FIG. 3: (a) Spectral response (in arbitrary units) for the dipolar excitation as a function of frequency uj (in units of uo), as 
calculated from the HF-RPA at T/Tc — 0.6 with the choice 77 = 0.005a;o. The density fluctuations at resonance (in arbitrary 
units) are plotted as functions of the radial coordinate r (in units of Uho = {h/muJoY^'^) in (b) for the condensate and in (c) for 
the non-condensate (solid lines). In the same panels are also shown the corresponding results from the analytical expressions 
of the mode eigenvectors (circles). 

thermal cloud near u) = 2.0a;o only a small change is present in the spectral intensity. 

V. DIPOLE MODE 

An important check on the accuracy of the theory is offered by the Kohn theorem. One can analytically prove 
that the dipolar oscillation in the a direction (with a = a;, y, or 2 in the general case of an anisotropic trap) 
is described by the Ansatz ~ {d/dva — mraUJa/fi) ^q, S^* = {d/dva + mraOJa/fi) Sh = dn^/dra, 5rh = 
{d/dva — 2mraOJa/K) fhP , and 5m* = {d/dva + 2mraUJa/h) fhP* . The theorem asserts that the corresponding mode 
frequency is given by the bare trap frequency uja- 

In Fig. 3(a) we show the spectral response for a dipolar excitation as obtained from the HF-RPA at T/Tc = 0.6 
and 77 = O.OOSa^Q. It has been explicitly shown that the Kohn theorem is satisfied in this approach 2(3,,^. As a result 
a sharp resonance is present in the HF-RPA dipole spectrum at lo = ujq. The density fluctuations at the resonance, 
as calculated from the solution of the dynamical equations, are plotted in Figs. 3(b) and 3(c) as solid lines and are 
compared with the predictions of the above Ansatz (circles). The two methods give almost the same result for both 
condensate and thermal density fluctuations, except for a weak structure in the thermal density fluctuation which 
may be due to the truncation of the basis sets. 

In Fig. 4 we show the spectral response of the dipole mode as obtained from the HFB-RPA with the same choice 
of parameters. In this approximation, the generalized Kohn theorem is not exactly satisfied, since a secondary peak 
is found in the spectrum at w ~ 1.13wo- According to the discussion given by Lewenstein and You p^ . a possible 
reason for this inaccuracy is the non-completeness of the set of quasiparticle wavefunctions used in the calculation. 
There also are appreciable distortions of the eigenvectors for the non-condensate oscillations in Fig. 4(c). 
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FIG. 4: The same as Fig. 3 for the HFB-RPA 
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FIG. 5; Spectral response (in arbitrary units) as a function of frequency lu (in units of ujo) for the monopole mode (a)and 
the quadrupole mode (b), as calculated with -q — O.OIljo from the HFB-RPA at the temperatures indicated in the figure. The 
curves are progressively shifted upwards by one unit for clarity and the quadrupole response at T/Tc = 0.1 is reduced by a 
factor of 3. The dashed line in each panel indicates how the condensate resonance moves with temperature. 



VI. MONOPOLE AND QUADRUPOLE MODES 



We present in this section the numerical results of the HFB-RPA for the monopole and quadrupole modes and 
compare them with those given by the SOBP-RPA and by the HF-RPA. These various theories give somewhat 
different results for the spectra at intermediate values of the temperature, in the range O.ATc <T< O.STc. 

In Fig. 5 we plot the HFB-RPA spectral functions at various temperatures. For ksT > two main resonances 
are seen in each spectrum, which can be interpreted as representing the collective oscillations of the non-condensate 
and of the condensate. The oscillator strength of each resonance has been extracted from the spectra and is shown in 
Fig. 6 as a function of temperature. Naturally, with increasing T/Tc the amplitude of the non-condensate resonances 
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FIG. 6: Amplitude of the HFB-RPA resonances (in arbitrary units) from Fig. 5 as a function of reduced temperature T/Tc 
for the monopole (a) and the quadrupole (b). The sohd and empty circles refer to the condensate and to the non-condensate, 
respectively. The lines are guides to the eye. 




CO/ COq CD/ CDq 

FIG. 7: Spectral response (in arbitrary units) as a function of frequency uj (in units of cjo) for the monopole mode (a) and the 
quadrupole mode (b), at T/Tc = 0.6 with rj = O.OlcJo from the HFB-RPA (solid lines), the SOBP-RPA (dashed lines) and the 
HF-RPA (dot-dashed lines). The arrows in panel (a) point to the condensate resonance position given by each RPA theory. 
The SOBP-RPA spectra as defined in Eqs. 11711 and 1181 do not include the contribution from the direct excitation of the 
non-condensate. 

grows (empty circles) while that of the condensate resonances decreases (solid circles). The amplitudes of the modes 
in the two components of the gas are comparable with each other near T/Tc — 0.5, where the non-condensate fraction 
is populated by about 30% for our choice of parameters. Above this temperature the strength of the non-condensate 
resonances increases very rapidly. 

In Fig. 7 we compare with each other the numerical results from the RPA theories for the monopolar and quadrupo- 
lar spectra at T/Tc — 0.6. We see that the HF-RPA and HFB-RPA closely agree in their predictions on the main 
non-condensate resonances for both types of excitations. We also see that all three theories predict essentially very 
similar results for the main quadrupolar resonance of the condensate, the position of the main peak at ~ 1.55a;o 
in Fig. 7(b) being also in agreement with the result of the HFB-Popov approximation (not shown). In the following 
we concentrate on the main condensate resonance in the monopolar mode, for which the three theories give rather 
different predictions as is emphasized by the three arrows in Fig. 7(a). In fact, the partial spectra of condensate 
and non-condensate show an appreciable overlap in this frequency range, implying a stronger dynamical coupling 
between the breathing excitations of the two components of the gas and therefore an enhanced sensitivity to the 
approximations made in the theory. 

To better illustrate the difference among the various theories, we extract the monopolar mode frequency of the 
condensate from the peak in x (^) ^.nd plot it in Fig. 8 as a function of reduced temperature. For comparison 
we also show the mode frequency given by the HFB-Popov theory (see Sec. II). The most remarkable feature of 
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FIG. 8: Monopole excitation frequency lom (in units of loq) as a function of reduced temperature T /Tc, as predicted by various 
theories: the HFB- Popov (solid line), the HFB-RPA (solid circles), the SOBP-RPA (empty circles) and the HF-RPA (stars). 
The lines connecting the symbols are guides to the eye. 




r r 



FIG. 9: Density fluctuations (in arbitrary units) as functions of the radial coordinate r (in units of aho) for the monopole mode 
(a) and the quadrupole mode (b), as calculated from the HFB-RPA for T /Tc — 0.6 at the appropriate excitation frequency of 
the condensate. In both panels the condensate density fluctuation is reduced by a factor of 5 for clarity. 



Fig. 8 is that all three RPA theories show a non-monotonic behavior of the resonance as a function of temperature, 
in contrast with the prediction of the HFB-Popov theory in which the resonance frequency decreases monotonically 
with increasing temperature. This difference is due to the dynamical coupling between the condensate and the non- 
condensate, which is neglected in the mean-field theory and becomes important as the non-condensate is significantly 
populated. 

Let us now compare the three RPA theories, which transcend the mean-field level. At low temperature (T/Tc < 0.4) 
we observe two different trends: the mode frequencies obtained from the HFB-RPA and from the SOBP-RPA are in 
close agreement and move upwards with temperature, whereas the mode frequency predicted by the HF-RPA tends 
to decrease. The latter trend is in good agreement with the HFB-Popov theory, in accord with the proof already 
given in Ref . . The upward trend of the mode frequency with temperature is manifested in all RPA theories at 
intermediate temperatures, reaching near T/Tc = 0.7 the highest sensitivity to the detailed description of the physical 
process in which the thermal cloud is driven by its self-generated dynamical potential. Finally, in proximity of the 
critical temperature all three theories tend to agree as the anomalous density fluctuations disappear. 
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The fact that a large upward frequency shift is found with increasing temperature in both the SOBP-RPA and the 
HFB-RPA suggests that a significant role is played by the anomalous density fluctuations. In Fig. 9 wc show the 
partial density fluctuations which accompany the monopolar and quadrupolar condensate resonances at T/Tc = 0.6, as 
calculated from the HFB-RPA. In both modes we find that the anomalous density fluctuations are at this temperature 
at least comparable in magnitude to the fluctuations of the normal density. 



VII. CONCLUSIONS 



In conclusion, we have developed a random-phase theory for the dynamics of a weakly interacting Bose gas under 
external confinement at finite temperature. In the theory the dynamics of the condensate and of the thermal cloud are 
treated on the same footing and a previous Hartree-Fock random-phase scheme is extended through the inclusion of 
the anomalous density fluctuations. The theory satisfies with good numerical accuracy the generalized Kohn theorem 
and correctly reduces to the second-order Beliaev-Popov theory if one neglects the process in which the thermal cloud 
is driven by its self-generated potential. It thereby fully includes the Landau-Beliaev damping mechanism. 

We have compared the theory with the second-order Beliaev-Popov theory and with the Hartree-Fock random- 
phase theory by numerical illustrations for a condensate of ^^Rb atoms inside a spherical trap. The location of the 
main monopolar and quadrupolar resonances of the thermal cloud are well reproduced in the Hartree-Fock RPA and 
the frequency of the quadrupole mode of the condensate does not differ significantly from the mean-field HFB-Popov 
prediction. We have instead found that for T > OATc the temperature dependence of the breathing mode frequency of 
the condensate obtained from the various RPA theories is very different from the HFB-Popov result. A significant role 
appears to be played in the dynamics of the Bosc-condensed gas by the anomalous density fluctuations of the thermal 
cloud at intermediate temperatures, even though they are known not to affect significantly the thermodynamics of 
the trapped gas [29ll30ll3l|. 

Our results, though restricted to isotropic confinement, may be relevant in connection with the JILA experiments 
, where the breathing mode in an anisotropic trap showed a frequency upshift with temperature which could not 
be accounted for by the HFB-Popov theory Q . A quantitative comparison between experimental data and the RPA 
predictions for an anisotropic trap would be interesting for a full test of the theory and we hope to address this issue 
in future work. 
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APPENDIX A: THE TWO-PARTICLE RESPONSE FUNCTIONS 



We present here a brief explanation on how to derive the response functions used in Eqs. © and (|10|) and list the 
two-particle response functions of the non-condensate. 

Let us consider for example the expression of Xcc (r,r';w). The most convenient way to obtain it is to calculate the 
bosonic Matsubara Green's function with imaginary time variable |22j . 

Xcc (r, r'; r) = - (Tr^{v, t)^{v' , 0))^ . (Al) 

Here T,- denote the ordering in imaginary time and (...)q denotes the equilibrium statistical average. By 
expressing the operator '0(r, r) in terms of the Bogoliubov quasiparticle operators cti and af ^ ^{r,T) = 
[uj{r)aje''^^'^ + v* {r)a'j e'^^'^] , we can rewrite Xcc (r, r'; r) in the form 

Xcc(r,r';T>0) = - (^(r, T)^+(r', 0))^ , 
j.k 

= K-(r)«;(r') (1 + /,) er^^- + v*{r)u,ir')f,e^^-] . (A2) 
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We then carry out a Fourier transform with respect to the imaginary time variable r, 







Xcc(r,r';io.„) = / dre'^-^Xcc (r, r'; r > 0) 



E 



u,{v)v*{v') v*{v)u,{v') 



(A3) 



where iujn ~ 2mri/(3. With the analytic continuation iujn — w + ijy we obtain the expression for Xcc (r,r';cLi) in Eq. 

CD- , 

The two-particle response functions of the non-condensate can be derived in a similar way. They take the following 
forms: 



with 



and 



with 



and 



with 



and 



Xnn (r, r'; uj) = (r, r'; uj) + x~fi (r, r'; uj) 



(A4) 



xis(i-,r';w) =^ 



{u*Uj + V*Vj) {u^U* + ViV*) if, - fj) 

huj'^ + (ei — ej) 



{U^VJ + V^Uj) {u*v* + v*u*) (1 + /, + fj) 



{u*V* + V*U*) {u^Vj + V^U■i) (1 + /, + fj) 

huj+ + (ei + f-j) 



Xnfa (r,r';w) = x^+s (r',!-;^) = xill (r, r'; w) xil (r,r';w) 



(A5) 



(1) . , ^ „ + ^Vj) UlV* ifi - fj) 

x'nJn (r,r;^) = 2}_^ ^ 



Xii(r,r';a;) = 2^ 



ViUjV*Vj {1 + ft + fj) U*V*UiUj {1 + fi + fj) 



hui^ — (ci -I- Cj) tiuj+ + {ei + ej 



Xhrh+ (r, r'; w) = x*rhn (r': ^) = xiVi+ l"^' ^) + Xsi,+ (i"' I"'; ^) 



(A6) 



xk^+(i",i-;w) = 2\ 



-I- (ei - ej) 



Xri.. (r,r';..) = 2^ 



ViUjU*U* (1 -I- /i + fj) U*V*jViVj (1 + /j + fj) 

fiLU+ - {ci + e-j) -t- (ei + ej) 



Xfr^rn (r, r'; cj) Xa+™+ (r', r; uj) = x^k (r, r'; uj) + x^k (r, r'; uj) 



(A7) 
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with 



(1) / / N , >r- (.A 



hu;+ + {ci - Cj) 



and 



ij 



UiUjV^v] {1 + fi + fj) V*V*UiUj {l + fi + fj)' 



TVjJ+ - [Ci + Cj 



Xmm+ (r, r'; w) = x!^k+ (•■' + X^k+ (''' ^ 



(A8) 



with 



X^^ik+ (r,r';a;)=4^. 



V*UjViU*j {fi - fj) 



and 



X^k+ (r,r';a;) = 2^ 



UiUjU*U* {1 + fi + fj) V*V*ViVj (1 + /i + /j) 



i '^j '-^J 



and finally 



Xm+m (r, r'; w) = x^+a (""^ ^''^ ^) + X^,Ia, ('^' "^'i ^) 



(A9) 



with 



(1) . / N , - /j) 

Xk+*(r,r;u;)=4}^. 



?ia;+ + (ci - Cj) 



and 



'^'i- (r,r';a;) =2y 



(1 + /i + /j) U>*Ui«j (1 + /, + /j) 



TujJ+ - (e, + Cj) 



7ia;+ + (e, + e^) 



In above expressions ti;+ = a; + iO+ and we have used abbreviations such as u*UjUiUj = u*{r)uj{r)ui{r')u*{r'), which 
means that in the product of four position-dependent functions the first two depend on r and the latter two on r'. 
X^^^ and Xo6^ in the above expressions correspond to the excitation of single thermal quasiparticles and of pairs of 
thermal quasiparticles, respectively. 
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